Začneme tým, že si do RStudia importujeme dáta. Ide o surové textové dáta, takže to nie je žiadny problém. Zvolíme import dát, zadáme cestu a prenesieme do nášho kódu:
hygdata_v3 <- read.csv("~/GitHub/HYG-Database/hygdata_v3.csv")
View(hygdata_v3)
Teraz už s nimi môžeme pracovať, najprv sa však pozrieme, ako sú usporiadané.
head(hygdata_v3)
Môžeme si všimnúť, že ide o pomerne komplexné dáta, ktoré sú bohaté okrem zaradenia do rôznych astronomických katalógov aj na údaje o relatívnej polohe a rýchlosti týchto vesmírnych objektov, ako aj o ich svietivosti a podobne. Máme teda na výber množstvo súvislostí, ktoré môžeme odhaľovať medzi jednotlivými parametrami.
Našou úplne prvou úlohou bude pozrieť sa, ako sú naše dáta v súbore rozdelené. Nebudeme sa trápiť s vykresľovaním jednotlivých charakteristík, radšej použijeme nástroje, ktoré sú v R už vbudované (summary). Zároveň by bolo ale fajn otestovať nedokonalosti v dátach, a ďalšiu informáciu nám môže poskytnúť aj modus, ktorý v bežnom summary nie je prítomný, a tak si ho pripravíme zvlášť a aplikujeme spolu so summary na naše dáta.
print(summary(hygdata_v3))
id hip hd hr gl bf
Min. : 0 Min. : 1 Min. : 1 Min. : 2 :115813 :116516
1st Qu.: 29903 1st Qu.: 29565 1st Qu.: 46734 1st Qu.:2281 GJ 1001 : 1 41The1Ori : 3
Median : 59807 Median : 59172 Median :110382 Median :4567 GJ 1002 : 1 100 Her: 2
Mean : 59807 Mean : 59170 Mean :114377 Mean :4564 GJ 1003 : 1 12Omi Cap : 2
3rd Qu.: 89710 3rd Qu.: 88763 3rd Qu.:175849 3rd Qu.:6850 GJ 1004 : 1 17Sig CrB : 2
Max. :119615 Max. :120404 Max. :358431 Max. :9110 GJ 1005A: 1 19 Lyn : 2
NA's :1658 NA's :20732 NA's :110589 (Other) : 3796 (Other) : 3087
proper ra dec dist pmra
:119468 Min. : 0.000 Min. :-89.782 Min. : 0.0 Min. :-4432.650
268 G. Cet: 1 1st Qu.: 6.217 1st Qu.:-36.421 1st Qu.: 115.1 1st Qu.: -15.460
3C 273 : 1 Median :12.128 Median : -1.640 Median : 213.7 Median : -1.680
82 G. Eri : 1 Mean :12.095 Mean : -1.985 Mean : 8764.8 Mean : -1.312
96 G. Psc : 1 3rd Qu.:18.117 3rd Qu.: 31.519 3rd Qu.: 390.6 3rd Qu.: 12.180
Acamar : 1 Max. :23.999 Max. : 89.569 Max. :100000.0 Max. : 6767.260
(Other) : 141
pmdec rv mag absmag spect
Min. :-5813.00 Min. :-386.9000 Min. :-26.70 Min. :-16.6800 K0 : 8569
1st Qu.: -22.39 1st Qu.: 0.0000 1st Qu.: 7.65 1st Qu.: 0.1390 G5 : 6021
Median : -5.76 Median : 0.0000 Median : 8.46 Median : 1.4950 A0 : 4943
Mean : -19.32 Mean : -0.2773 Mean : 8.43 Mean : 0.9921 F8 : 4366
3rd Qu.: 3.77 3rd Qu.: 0.0000 3rd Qu.: 9.17 3rd Qu.: 3.1590 G0 : 4237
Max. : 9999.99 Max. : 471.0000 Max. : 21.00 Max. : 19.6290 F5 : 3860
(Other):87618
ci x y z vx
Min. :-0.4000 Min. :-99950.39 Min. :-99979.25 Min. :-99964.98 Min. :-1.023e-01
1st Qu.: 0.3480 1st Qu.: -89.02 1st Qu.: -91.17 1st Qu.: -107.54 1st Qu.:-1.034e-05
Median : 0.6160 Median : -1.04 Median : -1.24 Median : -3.41 Median : 1.300e-07
Mean : 0.7115 Mean : -234.85 Mean : -39.57 Mean : -231.90 Mean :-2.891e-05
3rd Qu.: 1.0830 3rd Qu.: 86.27 3rd Qu.: 91.85 3rd Qu.: 94.97 3rd Qu.: 9.100e-06
Max. : 5.4600 Max. : 99982.37 Max. : 99996.07 Max. : 99862.51 Max. : 1.023e-01
NA's :1882
vy vz rarad decrad pmrarad
Min. :-1.023e-01 Min. :-1.023e-01 Min. :0.000 Min. :-1.56700 Min. :-2.149e-05
1st Qu.:-1.860e-06 1st Qu.:-1.998e-05 1st Qu.:1.628 1st Qu.:-0.63566 1st Qu.:-7.495e-08
Median : 1.182e-05 Median :-6.230e-06 Median :3.175 Median :-0.02862 Median :-8.140e-09
Mean : 2.164e-04 Mean :-1.566e-04 Mean :3.167 Mean :-0.03465 Mean :-6.360e-09
3rd Qu.: 3.209e-05 3rd Qu.: 3.150e-06 3rd Qu.:4.743 3rd Qu.: 0.55011 3rd Qu.: 5.905e-08
Max. : 1.023e-01 Max. : 1.023e-01 Max. :6.283 Max. : 1.56328 Max. : 3.281e-05
pmdecrad bayer flam con comp comp_primary
Min. :-2.818e-05 :118078 Min. : 1.0 Cen : 4222 Min. :1.000 Min. : 0
1st Qu.:-1.086e-07 Alp : 80 1st Qu.: 14.0 UMa : 3549 1st Qu.:1.000 1st Qu.: 29812
Median :-2.793e-08 Bet : 77 Median : 31.5 Her : 3383 Median :1.000 Median : 59628
Mean :-9.366e-08 Eps : 74 Mean : 37.3 Cyg : 3081 Mean :1.005 Mean : 59635
3rd Qu.: 1.828e-08 Del : 71 3rd Qu.: 55.0 Hya : 3027 3rd Qu.:1.000 3rd Qu.: 89453
Max. : 5.007e-05 Eta : 68 Max. :140.0 Cet : 2958 Max. :3.000 Max. :119615
(Other): 1166 NA's :116874 (Other):99394
base lum var var_min var_max
:118527 Min. : 0 :113624 Min. :-1.33 Min. :-1.52
GJ 1047 : 3 1st Qu.: 5 R : 60 1st Qu.: 8.53 1st Qu.: 8.24
Gl 57.1: 3 Median : 22 S : 52 Median : 9.85 Median : 9.65
Gl 60 : 3 Mean : 354469 T : 50 Mean : 9.50 Mean : 9.26
Gl 100 : 3 3rd Qu.: 77 W : 39 3rd Qu.:10.71 3rd Qu.:10.49
Gl 106.1: 3 Max. :409260660 U : 37 Max. :14.90 Max. :13.70
(Other) : 1072 (Other): 5752 NA's :102623 NA's :102623
Ak si pozorne prezrieme zobrazené charakteristiky, môžeme si všimnúť, že v niektorých položkách sú prítomné outlieri. Keďže však ide o reálne objekty s objektívnymi vlastnosťami, nebudeme ich zatiaľ vyradzovať (napríklad lum - svietivosť ako násobok svietivosti Slnka). Určite však budeme musieť upraviť veľmi dôležitú položku dist - teda vzdialenosť hviezdy získanú na základe jej paralaxy (relatívnej zmeny polohy kvôli pozorovaniu z dvoch protiľahlých miest na obežnej dráhe Slnka). Ako si môžeme všimnúť, tretí kvartil zasahuje hlboko pod maximálnu hodnotu súboru. Keď sa pozrieme na modus tohto spojitého súboru, zistíme, že je to práve táto maximálna hodnota. Táto maximálna presná hodnota v spojitom spektre je dosť podozrivá:
modus <- function(v)
{
uniqv <- unique(v)
uniqv[which.max(tabulate(match(v, uniqv)))]
}
modus(hygdata_v3[,10])
[1] 1e+05
hist(hygdata_v3[,10],xlab = "vzdialenosť(získaná pomocou paralaxy)",ylab = "počet výskytov",col = "darkblue",border = "blue")
Aha naozaj. Nájsť modus v spojitom súbore dát väčšinou znamená nejakú špecifickú vlastnosť značenia a histogram to potvrdzuje. V tomto bode si teda ešte raz prezrieme readme ku našej databáze. Naozaj. hodnota 100 000 znamená, že dáta chýbajú, alebo sú pochybné. Očistime naše dáta od tohoto znečistenia a pozrime sa na početnosti.
library(tidyverse)
[30m-- [1mAttaching packages[22m --------------------------------------- tidyverse 1.2.1 --[39m
[30m[32m<U+221A>[30m [34mggplot2[30m 3.2.1 [32m<U+221A>[30m [34mpurrr [30m 0.3.3
[32m<U+221A>[30m [34mtibble [30m 2.1.3 [32m<U+221A>[30m [34mdplyr [30m 0.8.3
[32m<U+221A>[30m [34mtidyr [30m 1.0.0 [32m<U+221A>[30m [34mstringr[30m 1.4.0
[32m<U+221A>[30m [34mreadr [30m 1.3.1 [32m<U+221A>[30m [34mforcats[30m 0.4.0[39m
[30m-- [1mConflicts[22m ------------------------------------------ tidyverse_conflicts() --
[31mx[30m [34mdplyr[30m::[32mfilter()[30m masks [34mstats[30m::filter()
[31mx[30m [34mdplyr[30m::[32mlag()[30m masks [34mstats[30m::lag()[39m
purified_1 <- filter(hygdata_v3,dist<=50000)
rest_1 <- filter(hygdata_v3, dist > 50000)
hist(purified_1[,10],xlab = "vzdialenosť(získaná pomocou paralaxy)",ylab = "počet výskytov",col = "darkblue",border = "blue")
hist(rest_1[,14],xlab = "magnitúda",ylab = "počet výskytov",col = "darkred",border = "gray")
Začnime s niečím jednoduchým. Radi by sme dokázali, že magnitúda (zdanlivá jasnosť hviezd) bude rásť (jej reprezentujúca číslica) so vzdialenosťou (prirodzený predpoklad v prípade homogénneho prostredia). Vynesme si teda tieto dve veličiny do grafu.
library(ggplot2)
ggplot(data = purified_1) +
geom_point(mapping = aes(x = dist, y = mag,size = 1/10^(mag),alpha = 1/100))
Z takéhoto grafu až tak veľa informácií nezískame. Dáta by bolo fajn si utriediť do nejakých skupín, a až takéto skupiny zobrazovať (a prípadne nastaviť vizuál podľa ich veľkostí). Poďme sa teda pustiť do našej prvej podúlohy - vyberme si z dát len zaujímavé položky a zoraďme ich podľa hodnoty.
hygdata_1 <- select(purified_1, id, proper, dist, rv, mag, absmag, con, lum)
hygdata_1
Teraz sa pokúsime nejako regulovať počet údajov. Jednotlivé objekty podľa vzdialenosti a magnitúdy zatriedime do vhodného počtu skupín. Využijeme zoradenie objektov vzostupne podľa vzdialenosťi, vykreslenie príslušnosti do skupín a funkciu cut.
dist
kat <- cut(hygdata_1$dist, 40, include.lowest=TRUE)
mutate(hygdata_1, kat)
hygdata_1
ggplot(data = hygdata_1) +
geom_point(mapping = aes(x = kat, y = mag, size = 1/10^(mag), alpha = 1/100), )
Ako si môžeme všimnúť, stále informačná hodnota takejto vizualizácie, nespĺňa naše požiadavky. Skúsime to inak: vizualizujme priamo závislosť vzdialenosti od magnitúry, ale príslušnosť ku jednotlivým vzialenostiam pre lepšiu prehľadnosť ofarbíme podľa rozdelenia spektra vzdialenosté do nastaviteľného počtu intervalov.
hygdata_1
kat <- cut(hygdata_1$dist, 25, include.lowest=TRUE, labels = FALSE)
mutate(hygdata_1, kat)
ggplot(data = hygdata_1) +
geom_point(mapping = aes(x = dist, y = mag, size = 1/10^(mag), alpha = 1/100), colour = kat)
Ako si môžeme všimnúť, vo funkcií cut sme použili labels = FALSE, pretože nechceme ručne zadávať potenciálne veľké množstvo názvov, ktoré nemajú informačnú hodnotu, len sú feed-om pre parameter colour pri kreslení grafu. To, že sme túto možnosť dali do kódu nám ale umožňuje voľne a nezáväzne adjustovať počet farebných hladín a počet “chlievikov”, ktoré zobrazujeme - vlastne to používame ako taký “kvázihistogram” - vidíme takto na našom grafe početnosti. Táto funkcia to delí na rovnaké podintervaly. Nám by pre vizualizáciu možno pomohlo, keby to delíme na kúsky s rovnakým počtom vzoriek - tak by sme získali informáciu o rozdelení našich dát a bolo by toto rozdelenie krásne vidno. Teraz sa to pokúsime spraviť. Hľadali sme možnosti až sme došli k jednoduchej variante príkazu cut, a to cut_number:
hygdata_1
kat <- cut_number(hygdata_1$dist, 20 , labels = FALSE)
mutate(hygdata_1, kat)
ggplot(data = hygdata_1) +
geom_point(mapping = aes(x = dist, y = mag, size = 1/1.1^(mag), alpha = 1/1000), colour = kat)
Teraz už sú dáta zobrazené pomerne prehľadne. Môžeme teda usúdiť, že náš predpoklad sa nepotvrdil. Vecou, ktorá je vďaka tejto vizualizácií najzrejmejšia, je rozdiel počtov hviezd v databáze pre jednotlivé vzdialenosti. Šírka stĺpcu (chlievika), ktorý obsahuje fixný počet hviezd sa so stúpajúcou vzdialenosťou zväčšuje - nie sme schopní detekovať všetky hviezdy (za predpokladu relatívnej homogenity rozmiestenenia hviezd vo vesmíre). Druhým viditeľným poznatkom je, že do vzdialenosti asi 50 parsekov sme schopní zachytiť aj objekty s veľmi malou zdanlivou magnitúdou. Potom ale počet takýchto objektov v databáze exponenciálne klesá. Tento pokles nie je spôsobený tým, že by vo vzdialenejších častiach menej jasné objekty neboli, ale tým, že ich nie sme schopní detekovať. V našej databáze teda tieto objekty nie sú. Teda ako to už býva, náš jednoduchý predpoklad “so stúpajúcou vzdialenosťou stúpa číslo magnitúdy” musíme korigovať ako “so stúpajúcou vzdialenosťou sme schopní zachytiť čoraz menej hviezd a žiadne s príliš veľkou magnitúdou (nízkou zdanlivou jasnosťou)”.
Druhým cieľom, ktorý by sme chceli splniť, je vykreslenie tvaru súhvezdí na oblohe podľa ich súradníc. Na takýto komplexnejšií problém využijeme fazety a identifikáciu podľa stĺpca con (skratka pre constellation). V prvom rade si z nášho komplexného dátového súboru vyberieme len položky, ktoré sú pre nás zaujímavé. Budeme postupovať analogicky ako v prvej časti.
hygdata_2 <- select(hygdata_v3, id, proper, mag, con, x, y, z, vx, vy, vz, rarad, decrad, pmrarad, pmdecrad)
hygdata_2 <- filter(hygdata_2,con!="")
Poďme sa pozrieť na rozdelenia jednotlivých súradníc - na základe tohoto sa rozhodneme, ktorú metódu vizualizácie zvolíme. Zároveň sa poďme pozrieť, koľko hodnôt môže nadobúdať kategorická premenná con - koľko je súhvezdí?
library(ggplot2)
ggplot(data.frame(hygdata_2$con), aes(x=hygdata_2$con), colour = hygdata_2$con ) +
geom_bar()
Ako môžeme vidieť, Súhvezdí je trochu moc, aby sme ich mohli všetky vykresliť. Vyberieme si teda len niektoré z nich a tejto vzorke sa budeme venovať. Samozrejme, v prípade záujmu by sme nasledujúce postupy mohli použiť aj na iných ľubovoľných súhvezdiach. Poďme sa pozrieť na počty hviezd v niektorých súhvezdiach / zvoľme napríklad Pisces - Ryby a Leo - lev.
hygdata_2Psc<-filter(hygdata_2, con =="Psc")
hist(hygdata_2Psc[,3],xlab = "magnitúda",ylab = "počet výskytov",col = "darkred",border = "gray")
hygdata_2Leo<-filter(hygdata_2, con=="Leo")
hist(hygdata_2Leo[,3],xlab = "magnitúda",ylab = "počet výskytov",col = "darkblue",border = "gray")
Ako vidíme, počty hviezd v jednotlivých súhvezdiach sú vyššie, ako by bolo možné prehľadne zobraziť. Poďme teda odfiltrovať hviezdy, ktoré nie sú viditeľné voľným okom. Uvádza sa, že voľným okom je možné vidieť objekty s veľkosťou okolo 6 magnitúdy (http://dust.ess.uci.edu/ppr/ppr_Wea47.pdf). Použime tento predpoklad a odfiltrujme hviezdy, ktoré nie sú viditeľné voľným okom (po optimalizácií sme zvolili ešte nižšiu hodnotu magnitúdy cutoff):
hygdata_2
hygdata_2_vi <- filter(hygdata_2, mag<4.5)
hygdata_2_vi
A môžeme začať vykresľovať! Samozrejme, ešte zvoľme, ktoré súhvezdia vykreslíme. Tak spravme túto úlohu pre 12 znamení zverokruhu:
hygdata_2_vi
hygdata_2_vi_z <- filter(hygdata_2_vi, con =="Ari"|con =="Tau"|con =="Gem"|con =="Cnc"|con =="Leo"|con =="Vir"|con =="Lib"|con =="Sco"|con =="Sgr"|con =="Aqr"|con =="Psc"|con =="Cap")
hygdata_2_vi_z
A teraz už naozaj môžeme začať vykresľovať. Skúsme, ako to vôbec bude vyzerať pre jedno súhvezdie (dobre identifikovateľný Leo). Vykresľujeme jednotlivé súhvezdia podľa ich súradníc na nebeskej klenbe (najprv sme experimentovali s x, y a z súradnicami až pokým sme si neuvedomili ich význam v tejto databáze). Takto sú charakteristické tvary viditeľné:
hygdata_2_vi_leo <- filter(hygdata_2_vi_z, con =="Leo")
ggplot(data = hygdata_2_vi_leo) +
geom_point(mapping = aes(x = hygdata_2_vi_leo$rarad, y = hygdata_2_vi_leo$decrad, size = 1/1.01^(hygdata_2_vi_leo$mag), alpha = 1))
Naozaj to funguje. Tvar súhvezdia, ktorý je bežne zobrazovaný, je zrkadlovo prevrátený (https://upload.wikimedia.org/wikipedia/commons/5/56/Leo_IAU.svg), ale v podstate ide o ten istý obrazec, teda sa nám podarilo vizualizovať jedno súhvezdie. Skúsme to spraviť pre všetky zvolené súhvezdia zo znamení zverokruhu. Ak volíme fazetové zobrazenie, je nutné nastaviť parameter scales ako voľný, pretože inak by sme nevideli nič (zlá mierka).
hygdata_2_vi_z
ggplot(data = hygdata_2_vi_z) +
geom_point(mapping = aes(x = hygdata_2_vi_z$rarad, y = hygdata_2_vi_z$decrad, size = 1/100^(hygdata_2_vi_z$mag), alpha = 1), colour = "black", background=2) +
facet_wrap(~con,scales = "free") + theme(panel.background = element_rect(fill = 'white', colour = 'white'))
Ignoring unknown parameters: background
Niektoré súhvezdia sú rozpoznateľné lepšie, iné horšie. Dalo by sa to zlepšiť napríklad iným nastavením filtračného parametra magnitúdy pre každé súhvezdie, podľa počtu hviezd. Niektoré súhvezdia totiž obsahujú aj hviezdy s menšou magnitúdou, ktoré sú ale význačné pre identifikáciu tvaru. My sa však ešte pozrieme na ďalšiu charakteristiku. Pozrieme sa, ako budú tieto súhvezdia vyzerať o 10 000 rokov. Takýto posun môžeme vizualizovať na základe radiálnych rýchlostí za rok, ktoré sú taktiež prítomné v databáze. Takže využime funkcie balíka dplyr:
hygdata_2_vi_z
hygdata_2_vi_z <- (mutate(hygdata_2_vi_z, rarad_10k = 100000*pmrarad + rarad, decrad_10k = 100000*pmdecrad + decrad ))
Teraz už zostáva tvar súhvezdí o 100 000 rokov len vykresliť. (nezabúdame, že je to pri uvažovanej konštantnej rýchlosti pohybu ako aj magnitúde, čo nie je pravidlom):
hygdata_2_vi_z
ggplot(data = hygdata_2_vi_z) +
geom_point(mapping = aes(x = hygdata_2_vi_z$rarad_10k, y = hygdata_2_vi_z$decrad_10k, size = 1/10^(hygdata_2_vi_z$mag), alpha = 1), colour = "black", background=2) +
facet_wrap(~con,scales = "free") + theme(panel.background = element_rect(fill = 'white', colour = 'white'))
Ignoring unknown parameters: background
Takto budú naše súhvezdia vyzerať o 100 000 rokov. Tvar niektorých zostane ešte rozpoznateľný, zatiaľ čo u iných bude deformovaný. (Pri parametri 1 000 000 rokov už sú súhvezdia nerozpoznateľné)
V tejto práci sme sa pokúsili využiť programovací jazyk R na organizáciu a vizualizáciu dát. Zamerali sme sa práve na tieto dve oblasti, pretože ide o najčastejšie používané typy práce s dátami. V prvom našom konkrétnom prípade sme formulovali hypotézu o vzťahu magnitúdy a vzdialenosti hviezdy, ktorú sme po preusporiadaní dát a ich vykreslení upravili. Dokázali sme tak, že vizualizácia a organizácia dát je veľmi dôležitá pre sformulovanie/potvrdenie/vyvrátenie vedeckej hypotézy. V druhej skôr zábavnej časti sme sa venovali vykresľovaniu zdanlivej polohy hviezd na oblohe podľa súhvezdí. Zredukovali sme výstupy z databázy, čo nám umožnilo rozumne vizualizovať dáta. Pridanou hodnotou bola prognóza, ktorá snáď budúcim astronómom uľahčí orientáciu na oblohe. :) Metódy by bolo možné zlepšiť použitím komplexnejších postupov.